/*
Eingangsvariablen/Zgr.:

        - Re ... Realteil des Signals
        - Im ... Imaginaerteil des Signals
        - sizeValues ... Anzahl der Abtastpunkte (:=sizeof(Re) bzw. sizeof(Im) bzw. sizeof(tsampled) )
        - Nfreq ... Größe des Frequenzarrays fdetect (:= sizeof(freq) ) 
        - freq ... das Frequenzarrays fdetect 
        - tsampled ... der Zeitvektor
        - start, ende ... Start-/Endpunkte der Integration (:=Summation im diskreten Fall)
        - P ... Zgr auf den vorab reservierten Speicher für das Spektrum  
        - Erg ... Ergebnisarray um den Index und den Wert der maximalen Frequenz zu speichern
        - sinArray ... Vorberechneten Sinuswerte
        - cosArray ... Vorberechneten Kosinuswerte

*******************************************************************************************************
Status: noch keine Vereinfachungen, blanke Implementierung.
*******************************************************************************************************/

#include"tS.h"

void NDFT(double* Re,double* Im,int sizeValues,double* freq, double* tsampled, int start, int ende, double* P_Re, double* P_Im,double* P_abs,double * Erg,double* sinArray, double* cosArray)
{
        double x,y,sumRe,sumIm;
        int i=0, k=0;
        Erg[0]=0;
        int tid;
        omp_set_num_threads(CORES);
        #pragma omp parallel private(x,y,sumRe,sumIm,tid,i) 
        {

        #ifdef CONTROL
                tid=omp_get_thread_num();
                #pragma omp critical    
                cout<<"Thread "<<tid<<" in NDFT sagt #HALLO#"<<endl;
        #endif

        #pragma omp for schedule(static) 
        for(k=start;k<=ende;k++)
        {
                sumRe=0;
                sumIm=0;
                for(i=0;i<sizeValues;i++)
                {
                        x = cosArray[(k-start)*sizeValues+i]/sqrt(sizeValues);
                        y = sinArray[(k-start)*sizeValues+i]/sqrt(sizeValues);
                        sumRe=sumRe+Re[i]*x+Im[i]*y;
                        sumIm=sumIm-Re[i]*y+x*Im[i];

                }
                P_Re[k-start]=sumRe;
                P_Im[k-start]=sumIm;
                P_abs[k-start]=sqrt(sumRe*sumRe+sumIm*sumIm);
                if(Erg[0] < P_abs[k-start])
                {
                        Erg[0]=P_abs[k-start];
                        Erg[1]=k-start;
                }
        }
        }                      //Warteschleife bis alle threads fertig sind......

}

//EOF

